The Data.

City-level information pulled from the following sources:

LA Almanac (race, age, and income) http://www.laalmanac.com/population/po38.php http://www.laalmanac.com/employment/em12.php

LA County vaccine percent by city http://publichealth.lacounty.gov/media/coronavirus/vaccine/vaccine-dashboard.htm

LA County age adjusted incidence rates http://dashboard.publichealth.lacounty.gov/covid19_surveillance_dashboard/

Data are incidence rates (proportions), vaccination rates (proportions), and characteristics by city.

Hypotheses:

All hypothesis are on the city-level and care must be taken when interpreting the results as they describe city-level variation and not individual level variation or co-variation.

Investigation of the December 2020-January 2021 peak for COVID 19 captured by a single day for incidence rates: 2021-01-07:

  1. How do race/ethnicity, age, and income describe variation in incidence rates?

Investigation of the August 2021 peak for COVID 19 captured by a single day for incidence rates: 2021-08-18:

  1. How do August vaccination proportions, race/ethnicity, age, income, and incidence rates for the Dec.-Jan. peak describe variation in incidence rates during the August peak for COVID 19?

Investigation of August 2021 vaccination rates for COVID 19 captured by a single day for vaccination rates: 2021-08-18:

  1. How do race/ethnicity, age, income, and incidence rates for the Dec.-Jan. peak describe variation in August vaccination rates?

Variable Distributions and Coding

Cities

d <- readRDS("data/analysis_data.rds") %>% 
  arrange(Date.x) %>% 
  filter(Date.x == "2021-08-18")

city.populations <- data.frame(d$city, d$population.x)
names(city.populations) <- c("City", "Population_Size")
city.populations <- city.populations[order(city.populations$Population_Size, decreasing=T), ]
kable(city.populations, format = "simple", align = 'c', caption = "Cities and Population Size")
Cities and Population Size
City Population_Size
121 Los Angeles 4043337
95 Santa Clarita 220424
44 Glendale 206493
62 Lancaster 161570
78 Palmdale 158968
82 Pomona 157869
108 Torrance 149268
35 East Los Angeles 125269
39 El Monte 117269
33 Downey 114263
52 Inglewood 113582
114 West Covina 108235
77 Norwalk 107622
18 Burbank 107180
25 Compton 99904
101 South Gate 98155
20 Carson 93846
97 Santa Monica 92446
48 Hawthorne 91301
119 Whittier 91216
4 Alhambra 86724
61 Lakewood 80362
15 Bellflower 77735
12 Baldwin Park 76769
68 Lynwood 72047
85 Redondo Beach 68697
11 Azusa 65963
26 Covina 65851
6 Arcadia 65735
42 Florence-Graham 64705
74 Montebello 64375
81 Pico Rivera 64284
75 Monterey Park 62262
43 Gardena 61310
51 Huntington Park 59484
104 South Whittier 59222
32 Diamond Bar 57515
80 Paramount 56023
46 Hacienda Heights 55926
88 Rosemead 55350
45 Glendora 52764
89 Rowland Heights 51022
22 Cerritos 50067
56 La Mirada 49599
5 Altadena 43620
14 Bell Gardens 43071
84 Rancho Palos Verdes 42747
73 Monrovia 42681
8 Westmont 42442
92 San Gabriel 40954
57 La Puente 40697
29 Culver City 39865
115 West Hollywood 36951
23 Claremont 36484
107 Temple City 36455
13 Bell 36332
70 Manhattan Beach 35999
58 La Verne 35322
120 Willowbrook 34913
16 Beverly Hills 34520
90 San Dimas 34516
63 Lawndale 33614
111 Walnut 30532
72 Maywood 28049
21 Castaic 27191
34 Duarte 26444
102 South Pasadena 26053
91 San Fernando 24612
28 Cudahy 24347
19 Calabasas 24323
76 East San Gabriel 24036
110 Valinda 23371
100 South El Monte 22680
64 Lennox 22542
113 West Carson 22086
105 Stevenson Ranch 20966
2 Agoura Hills 20883
67 Lomita 20729
54 La Crescenta-Montrose 19801
49 Hermosa Beach 19670
96 Santa Fe Springs 18364
7 Artesia 16795
40 El Segundo 16786
112 Walnut Park 16143
37 East Rancho Dominguez 15308
47 Hawaiian Gardens 14676
79 Palos Verdes Estates 13522
93 San Marino 13277
27 Charter Oak 13144
24 Commerce 13069
60 Lake Los Angeles 12994
69 Malibu 12961
83 Quartz Hill 12906
99 Signal Hill 11797
98 Sierra Madre 10989
116 West Puente Valley 9835
71 Marina del Rey 9411
103 South San Gabriel 8848
118 Westlake Village 8360
87 Rolling Hills Estates 8113
1 Acton 7971
59 Ladera Heights 7071
10 Avocado Heights 6775
36 East Pasadena 6403
106 Sun Village 6036
55 La Habra Heights 5455
38 East Whittier 5306
30 Del Aire 4393
3 Agua Dulce 4158
66 Littlerock 4021
9 Avalon 3869
109 Val Verde 3309
31 Desert View Highlands 2493
94 San Pasqual 2035
86 Rolling Hills 1940
50 Hidden Hills 1890
65 Leona Valley 1751
41 Elizabeth Lake 1661
53 Irwindale 1459
117 West Rancho Dominguez 1359
17 Bradbury 1069
City <- d$city

Incidence Rates

January Peak: 2021-01-07

Jan.IR <- d$adj_case_14day_rate.y
Jan.IR.cat <- relevel(as.factor(ifelse(Jan.IR > quantile(Jan.IR, probs=.75), "High",
                      ifelse(Jan.IR < quantile(Jan.IR, probs=.25), "Low",
                             "Med"))), ref="Med")

hist(Jan.IR, breaks=seq(from=0, to=round(max(Jan.IR),-3), by=200), col="red", main="January Incidence Rates for Cities", xlab="Incidence Rates")
abline(v=quantile(Jan.IR, probs=c(.25, .5, .75)), lwd=2, lty=2)

August Peak: 2021-08-18

Aug.IR <- d$adj_case_14day_rate.x

hist(Aug.IR, breaks=seq(from=0, to=round(max(Aug.IR),-2), by=100), col="red", main="August Incidence Rates for Cities", xlab="Incidence Rates")

Percent Vaccinated in Each City

Aug.Vax.percent <- d$dose1_all_c_prcent
Aug.Vax.cat <- relevel(as.factor(ifelse(Aug.Vax.percent > quantile(Aug.Vax.percent, probs=.75), "High",
                      ifelse(Aug.Vax.percent < quantile(Aug.Vax.percent, probs=.25), "Low",
                             "Med"))), ref="Med")

hist(Aug.Vax.percent, breaks=seq(from=0, to=100, by=5), col="blue", main="Percent Vaccination Rates for Cities", xlab="Percent Vaccinated")
abline(v=quantile(Aug.Vax.percent, probs=c(.25, .5, .75)), lwd=2, lty=2)

Race/Ethnicity

Proportion of Hispanic Individuals

CityHispanic <- d$Hispanic
CityHispanic.cat <- relevel(as.factor(ifelse(d$Hispanic > quantile(d$Hispanic, probs=.75), "High",
                      ifelse(d$Hispanic < quantile(d$Hispanic, probs=.25), "Low",
                             "Med"))), ref="Med")

hist(CityHispanic, breaks=seq(from=0, to=100, by=5), main="Proportion of Hispanic Individuals", xlab="Proportion")
abline(v=quantile(CityHispanic, probs=c(.25, .5, .75)), lwd=2, lty=2)

Proportion of Black Individuals

CityBlack <- d$Black
CityBlack.cat <- relevel(as.factor(ifelse(d$Black > 10, "High","Low")), "Low")


hist(CityBlack, breaks=seq(from=0, to=100, by=5), main="Proportion of Black Individuals", xlab="Proportion")
abline(v=10, lwd=2, lty=2)
text(10, 60, "Cutoff at 10%", pos=4)

Proportion of Asian Individuals

CityAsian <- d$Asian
CityAsian.cat <- relevel(as.factor(ifelse(d$Asian > quantile(d$Asian, probs=.75), "High",
                      ifelse(d$Asian < quantile(d$Asian, probs=.5), "Low",
                             "Med"))), ref="Low")

hist(CityAsian, breaks=seq(from=0, to=100, by=5), main="Proportion of Asian Individuals", xlab="Proportion")
abline(v=quantile(CityAsian, probs=c(.25, .5, .75)), lwd=2, lty=2)

Age

Age data is captured by the number/proportion of individuals within each city in a age category. Age categories include: “Age.Under.15”,“Age.15.17”,“Age.18.24”,“Age.25.34”,“Age.35.54”,“Age.55.64”,“Age.65.”

Three groups are created for the analysis:

  • Young: combined age categories for “Age.Under.15”,“Age.15.17”,“Age.18.24”

  • Middle: age category for “Age.35.54”

  • Old: combined age categories for “Age.55.64”,“Age.65.”

Age <- d[,c("Age.Under.15","Age.15.17","Age.18.24","Age.25.34","Age.35.54","Age.55.64","Age.65.")]
Age <- t(apply(Age, 1, FUN=function(v) { v/sum(v) }))

CityAgeYoung <- apply(Age[,c("Age.Under.15","Age.15.17","Age.18.24")], 1, sum)
CityAgeOld <- apply(Age[,c("Age.55.64","Age.65.")], 1, sum)
CityAgeMid <- Age[,"Age.35.54"]
CityAge <- cbind(CityAgeYoung, CityAgeOld)

CityAgeOld.cat <- relevel(as.factor(ifelse(CityAgeOld > quantile(CityAgeOld, probs=.75), "High","notHigh")), "notHigh")
CityAgeYoung.cat <- relevel(as.factor(ifelse(CityAgeYoung > quantile(CityAgeYoung, probs=.75), "High","notHigh")), "notHigh")

ftable(CityAgeOld.cat, CityAgeYoung.cat)
##                CityAgeYoung.cat notHigh High
## CityAgeOld.cat                              
## notHigh                              61   30
## High                                 30    0

Household Income

HouseholdIncome <- d$Households
HouseholdIncome.cat <- relevel(as.factor(ifelse(HouseholdIncome > quantile(HouseholdIncome, probs=.75), "High",
                      ifelse(HouseholdIncome < quantile(HouseholdIncome, probs=.25), "Low",
                             "Med"))), ref="Med")

hist(HouseholdIncome, breaks=seq(from=0, to=round(max(HouseholdIncome),-2), by=10000), main="Household Income Within Each City", xlab="Income", col="green")
abline(v=quantile(HouseholdIncome, probs=c(.25, .5, .75)), lwd=2, lty=2)

Variables for Analysis

Distributions

CityAge.m <- cbind(CityAgeOld.cat, CityAgeYoung.cat)
CityRaceEthnicty.m <- cbind(CityHispanic.cat, CityBlack.cat, CityAsian.cat)
d.analysis <- data.frame(Aug.IR, Jan.IR.cat, Aug.Vax.cat, Aug.Vax.percent, CityRaceEthnicty.m, CityAge.m, HouseholdIncome.cat)

summarytools::view(dfSummary(d.analysis, style = 'grid',
                               max.distinct.values = 10, plain.ascii = FALSE, valid.col = FALSE, headings = FALSE), method = "render")
No Variable Stats / Values Freqs (% of Valid) Graph Missing
1 Aug.IR [numeric] Mean (sd) : 436.5 (178.6) min < med < max: 39 < 447 < 1086 IQR (CV) : 186 (0.4) 110 distinct values 0 (0.0%)
2 Jan.IR.cat [factor] 1. Med 2. High 3. Low
61(50.4%)
30(24.8%)
30(24.8%)
0 (0.0%)
3 Aug.Vax.cat [factor] 1. Med 2. High 3. Low
61(50.4%)
30(24.8%)
30(24.8%)
0 (0.0%)
4 Aug.Vax.percent [numeric] Mean (sd) : 63.7 (11.2) min < med < max: 19.3 < 64.9 < 83.8 IQR (CV) : 14.6 (0.2) 104 distinct values 0 (0.0%)
5 CityHispanic.cat [integer] Mean (sd) : 1.7 (0.8) min < med < max: 1 < 1 < 3 IQR (CV) : 1 (0.5)
1:61(50.4%)
2:30(24.8%)
3:30(24.8%)
0 (0.0%)
6 CityBlack.cat [integer] Min : 1 Mean : 1.1 Max : 2
1:104(86.0%)
2:17(14.0%)
0 (0.0%)
7 CityAsian.cat [integer] Mean (sd) : 1.8 (0.8) min < med < max: 1 < 2 < 3 IQR (CV) : 2 (0.5)
1:60(49.6%)
2:30(24.8%)
3:31(25.6%)
0 (0.0%)
8 CityAgeOld.cat [integer] Min : 1 Mean : 1.2 Max : 2
1:91(75.2%)
2:30(24.8%)
0 (0.0%)
9 CityAgeYoung.cat [integer] Min : 1 Mean : 1.2 Max : 2
1:91(75.2%)
2:30(24.8%)
0 (0.0%)
10 HouseholdIncome.cat [factor] 1. Med 2. High 3. Low
61(50.4%)
30(24.8%)
30(24.8%)
0 (0.0%)

Generated by summarytools 0.9.9 (R version 4.1.0)
2021-11-04

Correlation

# need to make sure low is lowest when setting as.numeric (to make cor more clear)





d.num <- data.frame(Aug.IR, Jan.IR.cat, Aug.Vax.cat,
                    CityHispanic.cat, as.numeric(CityBlack.cat), as.numeric(CityAsian.cat), 
                    as.numeric(CityAgeOld.cat), as.numeric(CityAgeYoung.cat),
                    HouseholdIncome.cat) %>% 
  mutate(across(c(Jan.IR.cat, Aug.Vax.cat, CityHispanic.cat, HouseholdIncome.cat), 
                ~ as.numeric(relevel(.x, ref = "Low"))))

cormat <- cor(d.num, use="complete.obs", method="pearson")
corrplot(cormat, type="upper", order="original",
          col=brewer.pal(n=8, name="RdYlBu"),
          title = "",
          addCoef.col = "black",
          tl.cex=.5, number.cex=.5)

Pairwise Trend Plots

January IR vs. Proportion of Hispanic Individuals

CityHispanic.cont <- as.numeric(relevel(CityHispanic.cat, "Low"))
coef <- summary(lm(Jan.IR~CityHispanic.cont))$coef
ggplot(d.analysis, aes(x=CityHispanic.cont, y=Jan.IR)) +
  geom_point(color="blue", size=2) +
  geom_smooth(method=lm, color="black") +
  geom_text(x=2.5, y=3000, label=paste(expression(beta), "=", round(coef[2,1],2))) +
  geom_text(x=2.5, y=2500, label=paste("p-value", "=", round(coef[2,4],25))) +
  scale_x_continuous(name ="Proportion of Hispanic Individuals", breaks=c(1,2,3), labels=c("Low", "Med", "High"))+
  ylab("Jan Incidence Rates")
## `geom_smooth()` using formula 'y ~ x'

January IR vs. Household Income

coef <- summary(lm(Jan.IR~as.numeric(relevel(HouseholdIncome.cat, "Low"))))$coef
ggplot(d.analysis, aes(x=as.numeric(relevel(HouseholdIncome.cat, "Low")), y=Jan.IR)) +
  geom_point(color="blue", size=2) +
  geom_smooth(method=lm, color="black") +
  geom_text(x=2.5, y=3000, label=paste(expression(beta), "=", round(coef[2,1],2))) +
  geom_text(x=2.5, y=2500, label=paste("p-value", "=", round(coef[2,4],23))) +
  scale_x_continuous(name ="Household Income", breaks=c(1,2,3), labels=c("Low", "Med", "High"))+
  ylab("January Incidence Rates")
## `geom_smooth()` using formula 'y ~ x'

August IR vs. January IR

coef <- summary(lm(Aug.IR~Jan.IR))$coef
ggplot(d.analysis, aes(x=Jan.IR, y=Aug.IR)) +
  geom_point(color="blue", size=2) +
  geom_smooth(method=lm, color="black") +
  geom_text(x=3000, y=900, label=paste(expression(beta), "=", round(coef[2,1],2))) +
  geom_text(x=3000, y=800, label=paste("p-value", "=", round(coef[2,4],6)))+
  xlab("January Incidence Rates")+ylab("August Incidence Rates")
## `geom_smooth()` using formula 'y ~ x'

August IR vs. Proportion of Hispanic Individuals

CityHispanic.cont <- as.numeric(relevel(CityHispanic.cat, "Low"))
coef <- summary(lm(Aug.IR~CityHispanic.cont))$coef
ggplot(d.analysis, aes(x=CityHispanic.cont, y=Aug.IR)) +
  geom_point(color="blue", size=2) +
  geom_smooth(method=lm, color="black") +
  geom_text(x=2.5, y=900, label=paste(expression(beta), "=", round(coef[2,1],2))) +
  geom_text(x=2.5, y=800, label=paste("p-value", "=", round(coef[2,4],2))) +
  scale_x_continuous(name ="Proportion of Hispanic Individuals", breaks=c(1,2,3), labels=c("Low", "Med", "High"))+
  ylab("August Incidence Rates")
## `geom_smooth()` using formula 'y ~ x'

August IR vs. Household Income

coef <- summary(lm(Aug.IR~as.numeric(relevel(HouseholdIncome.cat, "Low"))))$coef
ggplot(d.analysis, aes(x=as.numeric(relevel(HouseholdIncome.cat, "Low")), y=Aug.IR)) +
  geom_point(color="blue", size=2) +
  geom_smooth(method=lm, color="black") +
  geom_text(x=2.5, y=900, label=paste(expression(beta), "=", round(coef[2,1],2))) +
  geom_text(x=2.5, y=800, label=paste("p-value", "=", round(coef[2,4],2))) +
  scale_x_continuous(name ="Household Income", breaks=c(1,2,3), labels=c("Low", "Med", "High"))+
  ylab("August Incidence Rates")
## `geom_smooth()` using formula 'y ~ x'

August IR vs. August Percent Vaccinated

coef <- summary(lm(Aug.IR~Aug.Vax.percent))$coef
ggplot(d.analysis, aes(x=Aug.Vax.percent, y=Aug.IR)) +
  geom_point(color="blue", size=2) +
  geom_smooth(method=lm, color="black") +
  geom_text(x=80, y=900, label=paste(expression(beta), "=", round(coef[2,1],2))) +
  geom_text(x=80, y=800, label=paste("p-value", "=", round(coef[2,4],2)))+
  xlab("August Vaccination Percent")+ylab("August Incidence Rates")
## `geom_smooth()` using formula 'y ~ x'

January Incident Rates:

Regression Model

reg.Jan.noVax <- lm(Jan.IR ~ HouseholdIncome.cat+CityRaceEthnicty.m+CityAge.m)

Regression summary

kable(summary(reg.Jan.noVax)$coef, digits=3, format = "simple", align = 'c', caption = "January Regression")
January Regression
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1087.329 431.053 2.522 0.013
HouseholdIncome.catHigh -983.811 171.194 -5.747 0.000
HouseholdIncome.catLow 783.304 163.598 4.788 0.000
CityRaceEthnicty.mCityHispanic.cat 97.188 81.787 1.188 0.237
CityRaceEthnicty.mCityBlack.cat 65.163 165.618 0.393 0.695
CityRaceEthnicty.mCityAsian.cat 67.657 73.121 0.925 0.357
CityAge.mCityAgeOld.cat -107.726 146.655 -0.735 0.464
CityAge.mCityAgeYoung.cat 444.903 164.853 2.699 0.008
#stargazer(reg.Jan.noVax, type = 'html', ci=TRUE, ci.level=0.95, single.row = T)

Regression ANOVA

\(R^2 =\) 0.62

kable(anova(reg.Jan.noVax), digits=3, format = "simple", align = 'c', caption = "January Regression ANOVA")
January Regression ANOVA
Df Sum Sq Mean Sq F value Pr(>F)
HouseholdIncome.cat 2 60973275.2 30486637.6 87.403 0.000
CityRaceEthnicty.m 3 542565.2 180855.1 0.518 0.670
CityAge.m 2 3052626.7 1526313.3 4.376 0.015
Residuals 113 39415190.7 348807.0 NA NA

August Incident Rates:

Vaccination Only

Regression Model

reg.Aug.vax <- lm(Aug.IR ~ Aug.Vax.cat)

Regression summary

kable(summary(reg.Aug.vax)$coef, digits=3, format = "simple", align = 'c', caption = "August Regression")
August Regression
Estimate Std. Error t value Pr(>|t|)
(Intercept) 450.118 21.753 20.692 0.000
Aug.Vax.catHigh -109.075 37.886 -2.879 0.005
Aug.Vax.catLow 54.182 37.886 1.430 0.155

Regression ANOVA

\(R^2 =\) 0.11

kable(anova(reg.Aug.vax), digits=3, format = "simple", align = 'c', caption = "August Regression ANOVA")
August Regression ANOVA
Df Sum Sq Mean Sq F value Pr(>F)
Aug.Vax.cat 2 422581 211290.5 7.32 0.001
Residuals 118 3406082 28865.1 NA NA

Full Model

Regression Model

reg.Aug <- lm(Aug.IR ~ Aug.Vax.cat + Jan.IR.cat + HouseholdIncome.cat+CityRaceEthnicty.m+CityAge.m)

Regression summary

kable(summary(reg.Aug)$coef, digits=3, format = "simple", align = 'c', caption = "August Regression")
August Regression
Estimate Std. Error t value Pr(>|t|)
(Intercept) 298.772 114.632 2.606 0.010
Aug.Vax.catHigh -53.373 42.644 -1.252 0.213
Aug.Vax.catLow 43.507 41.916 1.038 0.302
Jan.IR.catHigh -37.110 47.008 -0.789 0.432
Jan.IR.catLow -218.582 49.835 -4.386 0.000
HouseholdIncome.catHigh 59.182 48.820 1.212 0.228
HouseholdIncome.catLow -67.747 48.738 -1.390 0.167
CityRaceEthnicty.mCityHispanic.cat 73.615 25.706 2.864 0.005
CityRaceEthnicty.mCityBlack.cat 88.676 44.028 2.014 0.046
CityRaceEthnicty.mCityAsian.cat -5.201 20.844 -0.250 0.803
CityAge.mCityAgeOld.cat -70.434 39.484 -1.784 0.077
CityAge.mCityAgeYoung.cat 58.694 45.717 1.284 0.202

Regression ANOVA

\(R^2 =\) 0.31

kable(anova(reg.Aug), digits=3, format = "simple", align = 'c', caption = "August Regression ANOVA")
August Regression ANOVA
Df Sum Sq Mean Sq F value Pr(>F)
Aug.Vax.cat 2 422581.0 211290.50 8.714 0.000
Jan.IR.cat 2 233137.6 116568.79 4.808 0.010
HouseholdIncome.cat 2 98384.8 49192.40 2.029 0.136
CityRaceEthnicty.m 3 287587.3 95862.42 3.954 0.010
CityAge.m 2 144166.4 72083.21 2.973 0.055
Residuals 109 2642806.1 24245.93 NA NA

August Vaccination Rates:

Regression Model with August Vaccination

reg.Aug.vax <- lm(Aug.Vax.percent ~ Jan.IR.cat+HouseholdIncome.cat+CityRaceEthnicty.m+CityAge.m)

Regression summary

kable(summary(reg.Aug.vax)$coef, digits=3, format = "simple", align = 'c', caption = "August Vaccination Regression")
August Vaccination Regression
Estimate Std. Error t value Pr(>|t|)
(Intercept) 59.663 6.915 8.627 0.000
Jan.IR.catHigh -1.763 2.791 -0.632 0.529
Jan.IR.catLow 0.389 2.866 0.136 0.892
HouseholdIncome.catHigh 3.236 2.958 1.094 0.276
HouseholdIncome.catLow -0.898 2.940 -0.305 0.761
CityRaceEthnicty.mCityHispanic.cat 3.718 1.521 2.445 0.016
CityRaceEthnicty.mCityBlack.cat 0.993 2.651 0.375 0.709
CityRaceEthnicty.mCityAsian.cat 3.690 1.176 3.137 0.002
CityAge.mCityAgeOld.cat -2.160 2.381 -0.907 0.366
CityAge.mCityAgeYoung.cat -6.132 2.681 -2.287 0.024

Regression ANOVA

kable(anova(reg.Aug.vax), digits=3, format = "simple", align = 'c', caption = "August Vaccination Regression ANOVA")
August Vaccination Regression ANOVA
Df Sum Sq Mean Sq F value Pr(>F)
Jan.IR.cat 2 1976.111 988.056 11.094 0.000
HouseholdIncome.cat 2 1059.438 529.719 5.948 0.004
CityRaceEthnicty.m 3 1696.679 565.560 6.350 0.001
CityAge.m 2 488.550 244.275 2.743 0.069
Residuals 111 9885.456 89.058 NA NA

Multnomial Regression Model with August Vaccination

reg.Aug.vax <- lm(Aug.Vax.percent ~ Jan.IR.cat+HouseholdIncome.cat+CityRaceEthnicty.m+CityAge.m)
reg.Aug.vax <- multinom(Aug.Vax.cat ~ Jan.IR.cat+HouseholdIncome.cat+CityRaceEthnicty.m+CityAge.m)
## # weights:  33 (20 variable)
## initial  value 132.932087 
## iter  10 value 78.668921
## iter  20 value 76.286634
## iter  30 value 76.182929
## iter  40 value 76.181907
## final  value 76.181903 
## converged

Regression summary for High Vaccination Rates

coef.high <- summary(reg.Aug.vax)$coefficients["High",]
OR.high <- exp(coef.high)
se.high <- summary(reg.Aug.vax)$standard.errors["High",]
z.high <- coef.high/se.high
p.value.high <- (1 - pnorm(abs(z.high), 0, 1)) * 2
r.high <- data.frame(OR.high, coef.high, se.high, z.high, p.value.high)
kable(r.high, digits=3, format = "simple", align = 'c', caption = "August Vaccination Regression for High Vax Rates")
August Vaccination Regression for High Vax Rates
OR.high coef.high se.high z.high p.value.high
(Intercept) 459743.494 13.038 1.130 11.538 0.000
Jan.IR.catHigh 0.000 -16.605 0.000 -12741324.270 0.000
Jan.IR.catLow 6.620 1.890 0.742 2.548 0.011
HouseholdIncome.catHigh 0.783 -0.245 0.834 -0.294 0.769
HouseholdIncome.catLow 2.879 1.058 1.343 0.787 0.431
CityRaceEthnicty.mCityHispanic.cat 1.589 0.463 0.413 1.122 0.262
CityRaceEthnicty.mCityBlack.cat 0.000 -15.583 1.130 -13.790 0.000
CityRaceEthnicty.mCityAsian.cat 1.028 0.028 0.392 0.071 0.944
CityAge.mCityAgeOld.cat 1.604 0.473 0.682 0.693 0.488
CityAge.mCityAgeYoung.cat 1.068 0.065 1.265 0.052 0.959

Regression summary for Low Vaccination Rates

coef.low <- summary(reg.Aug.vax)$coefficients["Low",]
OR.low <- exp(coef.low)
se.low <- summary(reg.Aug.vax)$standard.errors["Low",]
z.low <- coef.low/se.low
p.value.low <- (1 - pnorm(abs(z.low), 0, 1)) * 2
r.low <- data.frame(OR.low, coef.low, se.low, z.low, p.value.low)
kable(r.low, digits=3, format = "simple", align = 'c', caption = "August Vaccination Regression for Low Vax Rates")
August Vaccination Regression for Low Vax Rates
OR.low coef.low se.low z.low p.value.low
(Intercept) 3.841 1.346 2.349 0.573 0.567
Jan.IR.catHigh 0.791 -0.234 0.942 -0.249 0.803
Jan.IR.catLow 4.193 1.433 1.110 1.291 0.197
HouseholdIncome.catHigh 0.533 -0.629 1.031 -0.610 0.542
HouseholdIncome.catLow 1.750 0.560 0.934 0.600 0.549
CityRaceEthnicty.mCityHispanic.cat 0.377 -0.976 0.574 -1.701 0.089
CityRaceEthnicty.mCityBlack.cat 0.861 -0.149 0.708 -0.211 0.833
CityRaceEthnicty.mCityAsian.cat 0.081 -2.516 0.799 -3.151 0.002
CityAge.mCityAgeOld.cat 2.407 0.878 0.853 1.029 0.303
CityAge.mCityAgeYoung.cat 3.720 1.314 0.807 1.629 0.103




Model params over time

Full model fit by date
{r} adj_case_14day_rate.x ~ Vax.percent + Jan.IR.cat + HouseholdIncome.cat+CityHispanic.cat+CityBlack.cat+ CityAsian.cat+ CityAgeOld.cat+ CityAgeYoung.cat
No vaccination lag time in this example

Vaccination rates

Modeled as a continuous variable. Plots include dates > 04/01/2021

d <- readRDS("data/analysis_data.rds") 
lac_summarytable_dash <- read.csv("data/LA_County_Covid19_CSA_14day_case_death_table_Clean.csv")

# data clean ----
dout <- d %>% 
  mutate(Date.x = ymd(Date.x), 
         Vax.percent = dose1_all_c_prcent, 
         Jan.IR =  as.numeric(adj_case_14day_rate.y),
         Jan.IR.cat =  relevel(as.factor(ifelse(Jan.IR > quantile(Jan.IR, probs=.75), "High",
                                                ifelse(Jan.IR < quantile(Jan.IR, probs=.25), "Low",
                                                       "Med"))), ref="Low")) %>%
  dplyr::select(Date.x, city, adj_case_14day_rate.x, 
                Vax.percent, Jan.IR.cat, 
                HouseholdIncome.cat, CityHispanic.cat, CityBlack.cat, 
                CityAsian.cat, CityAgeOld.cat, CityAgeYoung.cat)




test <- dout %>% 
  # dplyr::filter(Date.x > "2021-03-01") %>% 
  tidyr::nest(data = -Date.x) %>% 
  dplyr::mutate(fit = purrr::map(data, ~lm(adj_case_14day_rate.x ~ Vax.percent + Jan.IR.cat + HouseholdIncome.cat+CityHispanic.cat+CityBlack.cat+ CityAsian.cat+ CityAgeOld.cat+ CityAgeYoung.cat, data = .x)),
                tidied = purrr::map(fit, ~broom::tidy(.x), quality = purrr::map(fit, ~glance(.x)))) %>% 
  dplyr::select(-data, -fit) %>% 
  tidyr::unnest(tidied) %>% 
  # dplyr::filter(grepl("Vax", term)) %>%
  arrange(Date.x)



# Vaccination rates ---- 
# after april 2021 (vaccination rates too low prior)
p1 <- ggplot(test %>% filter(grepl("Vax.percent", term), Date.x > "2021-04-01"), aes(Date.x, estimate)) + 
  geom_point() + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error)) +
  geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))
  # geom_line(data = lac_summarytable_dash %>% filter(city == "Los Angeles County"), aes(Date, adj_case_14day_rate))

p2 <- ggplot(lac_summarytable_dash %>% 
               filter(city == "Los Angeles County", 
                      Date > "2021-04-01") %>% 
               mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
  geom_line()

plot_grid(p1, p2, ncol = 1, align = 'v')

Race

p1 <- ggplot(test %>% filter(term %in% c("CityAsian.catHigh", "CityBlack.catHigh", "CityHispanic.catHigh")), aes(Date.x, estimate)) + 
  geom_point(aes(colour = term)) + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
  geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))

p2 <- ggplot(lac_summarytable_dash %>% 
               filter(city == "Los Angeles County") %>% 
               mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
  geom_line()

plot_grid(p1, p2, ncol = 1, align = 'v')

Household income

p1 <- ggplot(test %>% filter(term %in% c("HouseholdIncome.catHigh", "HouseholdIncome.catLow")), aes(Date.x, estimate)) + 
  geom_point(aes(colour = term)) + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
  geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.2))

p2 <- ggplot(lac_summarytable_dash %>% 
               filter(city == "Los Angeles County") %>% 
               mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
  geom_line()

plot_grid(p1, p2, ncol = 1, align = 'v')
## Warning: Removed 1 row(s) containing missing values (geom_path).

Age category

p1 <- ggplot(test %>% filter(term %in% c("CityAgeOld.catHigh", "CityAgeYoung.catHigh")), aes(Date.x, estimate)) + 
  geom_point(aes(colour = term)) + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
  geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))

p2 <- ggplot(lac_summarytable_dash %>% 
               filter(city == "Los Angeles County") %>% 
               mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
  geom_line()

plot_grid(p1, p2, ncol = 1, align = 'v')
## Warning: Removed 1 row(s) containing missing values (geom_path).

January incidence rates

p1 <- ggplot(test %>% filter(term %in% c("Jan.IR.catHigh", "Jan.IR.catMed")), aes(Date.x, estimate)) + 
  geom_point(aes(colour = term)) + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
  geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))
# geom_line(data = lac_summarytable_dash %>% filter(city == "Los Angeles County"), aes(Date, adj_case_14day_rate))
# ggsave(p1, width = 10, height = 8, file = "~/Desktop/test1.png")
p2 <- ggplot(lac_summarytable_dash %>% 
               filter(city == "Los Angeles County") %>% 
               mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
  geom_line()

plot_grid(p1, p2, ncol = 1, align = 'v')
## Warning: Removed 1 row(s) containing missing values (geom_path).




Model params over time (vax lag)

Full model fit by date
{r} adj_case_14day_rate.x ~ Vax.percent + Jan.IR.cat + HouseholdIncome.cat+CityHispanic.cat+CityBlack.cat+ CityAsian.cat+ CityAgeOld.cat+ CityAgeYoung.cat
Introduce 14 day lag between vaccination rates and COVID incidence rates

Vaccination rates

Modeled as a continuous variable. Plots include dates > 04/01/2021

d <- readRDS("data/analysis_data_14vax_lag.rds") 
lac_summarytable_dash <- read.csv("data/LA_County_Covid19_CSA_14day_case_death_table_Clean.csv")

# data clean ----
dout <- d %>% 
  mutate(Date.x = ymd(Date.x), 
         Vax.percent = dose1_all_c_prcent, 
         Jan.IR =  as.numeric(adj_case_14day_rate.y),
         Jan.IR.cat =  relevel(as.factor(ifelse(Jan.IR > quantile(Jan.IR, probs=.75), "High",
                                                ifelse(Jan.IR < quantile(Jan.IR, probs=.25), "Low",
                                                       "Med"))), ref="Low")) %>%
  dplyr::select(Date.x, city, adj_case_14day_rate.x, 
                Vax.percent, Jan.IR.cat, 
                HouseholdIncome.cat, CityHispanic.cat, CityBlack.cat, 
                CityAsian.cat, CityAgeOld.cat, CityAgeYoung.cat)




test <- dout %>% 
  # dplyr::filter(Date.x > "2021-03-01") %>% 
  tidyr::nest(data = -Date.x) %>% 
  dplyr::mutate(fit = purrr::map(data, ~lm(adj_case_14day_rate.x ~ Vax.percent + Jan.IR.cat + HouseholdIncome.cat+CityHispanic.cat+CityBlack.cat+ CityAsian.cat+ CityAgeOld.cat+ CityAgeYoung.cat, data = .x)),
                tidied = purrr::map(fit, ~broom::tidy(.x), quality = purrr::map(fit, ~glance(.x)))) %>% 
  dplyr::select(-data, -fit) %>% 
  tidyr::unnest(tidied) %>% 
  # dplyr::filter(grepl("Vax", term)) %>%
  arrange(Date.x)



# Vaccination rates ---- 
# after april 2021 (vaccination rates too low prior)
p1 <- ggplot(test %>% filter(grepl("Vax.percent", term), Date.x > "2021-04-01"), aes(Date.x, estimate)) + 
  geom_point() + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error)) +
  geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))
  # geom_line(data = lac_summarytable_dash %>% filter(city == "Los Angeles County"), aes(Date, adj_case_14day_rate))

p2 <- ggplot(lac_summarytable_dash %>% 
               filter(city == "Los Angeles County", 
                      Date > "2021-04-01") %>% 
               mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
  geom_line()

plot_grid(p1, p2, ncol = 1, align = 'v')

Race

p1 <- ggplot(test %>% filter(term %in% c("CityAsian.catHigh", "CityBlack.catHigh", "CityHispanic.catHigh")), aes(Date.x, estimate)) + 
  geom_point(aes(colour = term)) + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
  geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))

p2 <- ggplot(lac_summarytable_dash %>% 
               filter(city == "Los Angeles County") %>% 
               mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
  geom_line()

plot_grid(p1, p2, ncol = 1, align = 'v')

Household income

p1 <- ggplot(test %>% filter(term %in% c("HouseholdIncome.catHigh", "HouseholdIncome.catLow")), aes(Date.x, estimate)) + 
  geom_point(aes(colour = term)) + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
  geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.2))

p2 <- ggplot(lac_summarytable_dash %>% 
               filter(city == "Los Angeles County") %>% 
               mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
  geom_line()

plot_grid(p1, p2, ncol = 1, align = 'v')
## Warning: Removed 1 row(s) containing missing values (geom_path).

Age category

p1 <- ggplot(test %>% filter(term %in% c("CityAgeOld.catHigh", "CityAgeYoung.catHigh")), aes(Date.x, estimate)) + 
  geom_point(aes(colour = term)) + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
  geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))

p2 <- ggplot(lac_summarytable_dash %>% 
               filter(city == "Los Angeles County") %>% 
               mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
  geom_line()

plot_grid(p1, p2, ncol = 1, align = 'v')
## Warning: Removed 1 row(s) containing missing values (geom_path).

January incidence rates

p1 <- ggplot(test %>% filter(term %in% c("Jan.IR.catHigh", "Jan.IR.catMed")), aes(Date.x, estimate)) + 
  geom_point(aes(colour = term)) + 
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error, colour = term)) +
  geom_hline(yintercept = 0) + theme_bw() + theme(legend.position = c(0.8, 0.8))
# geom_line(data = lac_summarytable_dash %>% filter(city == "Los Angeles County"), aes(Date, adj_case_14day_rate))
# ggsave(p1, width = 10, height = 8, file = "~/Desktop/test1.png")
p2 <- ggplot(lac_summarytable_dash %>% 
               filter(city == "Los Angeles County") %>% 
               mutate(Date = ymd(Date)), aes(Date, adj_case_14day_rate)) +
  geom_line()

plot_grid(p1, p2, ncol = 1, align = 'v')
## Warning: Removed 1 row(s) containing missing values (geom_path).